function [idat outp] = fsb_realign_dft3D(idat,param)

% FSB - DEV: Realign scans with discrete fourier transform algorithm volumewise
%
% EXAMPLE:
% [idat outp] = fsb_realign_dft3D(idat,1)
%
% INPUT:
% idat:         4-D image data
% ref:          Reference slice
%
% OUTPUT:
% A 4-D image data volume, realigned in the y direction only
%
% CALLED BY:
% fsb_realign.m
%
% NOTES:
% Requires the dftregistration function copyright
% Manuel Guizar available on matlabcentral.com
% (efficient subpixel registration)
% Citation for his algorithm:
% Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
% "Efficient subpixel image registration algorithms," Opt. Lett. 33,
% 156-158 (2008).
% This function uses the same algorithm as the fsb_realign_dft function,
% but does not realign slice by slice. It rather calculates the mean of the
% slice shift values and then corrects for this in a second step
% Precision is coarser than latter function.
% Averaging between subsequent volumes creates additional smoothness
% This is due to the algorithm used for shifting images back into position
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if second input is provided
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if nargin<2
    param = 0.01;
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Go over volumes and slices and check for mean movement in y direction
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

reafig = waitbar(0,'Realigning ...');
[a b c d] = size(idat);
midat = single(mean(idat,4));
tic;
for y=1:c;
    idat1 = single(squeeze(midat(:,:,y)));
    waitbar (y/c);
    for x = 1:d;
        idatx = single(squeeze(idat(:,:,y,x)));
        output = dftregistration(fft2(idat1),fft2(idatx),param);
        outp.output(x,:,y) = output;
    end
end

shiftpar = double(squeeze(outp.output(:,1,:)));
meanshift = mean(shiftpar,2);
dca = fix(meanshift);
disp (['Max. negative translation: ' num2str(min(meanshift)) ' voxel']);
disp (['Max. positive translation: ' num2str(max(meanshift)) ' voxel']);
dc_rem = meanshift-dca;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Shift every single volume back or forth by detected movement
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Shifting volumes...');
for x = 1:d

    %~~~~~~~~~~~~~~~~~~~~
    % Create dummy volume
    %~~~~~~~~~~~~~~~~~~~~
    p_idat = idat(:,:,:,x);

    %~~~~~~~~~~~~~~~~~~~~
    % Shift volumes if necessary
    %~~~~~~~~~~~~~~~~~~~~
    if dca(x)~=0 && ~isnan(dca(x));
        p_idat = circshift(p_idat,[0 (-dca(x))]);
    end

    %~~~~~~~~~~~~~~~~~~~~
    % Create shifted volumes for averaging
    %~~~~~~~~~~~~~~~~~~~~
    p_idat1 = circshift(p_idat,[0 1]);
    p_idat2 = circshift (p_idat,[0 -1]);

    %~~~~~~~~~~~~~~~~~~~~
    % Average original and shifted volume to approximate realignment
    %~~~~~~~~~~~~~~~~~~~~
    if dc_rem(x)>0;
        p_idat = (p_idat2*dc_rem(x))+(p_idat*(1-dc_rem(x)));
    elseif dc_rem(x)<0;
        p_idat = (p_idat1*abs(dc_rem(x)))+(p_idat*(1-abs(dc_rem(x))));
    end

    %~~~~~~~~~~~~~~~~~~~~
    % write into idat
    %~~~~~~~~~~~~~~~~~~~~
    idat(:,:,:,x) = p_idat;

    waitbar(x/d);
end
close (h);

toc;
close (reafig);
end


